和part1不同的是control不一样
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]
file2=file2[!is.na(file2$unitID),]

con=read.table("./20201201/all.site",head=T,sep="\t")
con$id=paste(con$Chr,con$Pos,sep=":")
con=con[!con$id %in% as.character(filea$id),]		###remove the ASH

asm=read.table("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS/ASM_log/220520ASM_bias_dscore_all.txt",head=TRUE,sep="\t")
asm$unitID=paste(asm$chr,asm$start,sep=":")
length(unique(asm$unitID))
asmid=unique(as.character(asm$unitID))

###AJHG
id1=as.character(file2$id)
ajhg=read.table("E:/0 公共数据库差异情况/db_for_5hmc/2704_ASMs_hg18ToHg19_AJHG.bed",head=F,sep="\t")
ajhg$id=paste(ajhg$V1,ajhg$V2,sep=":")

length(intersect(as.character(file$id),as.character(ajhg$id)))

length(intersect(as.character(filea$id),as.character(asmid)))

c=length(intersect(as.character(con$id),as.character(ajhg$id)))
d=dim(con)[1]-c

result=data.frame(matrix(,nrow = 5,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",c,NA,NA)

a=length(intersect(as.character(filea$id),as.character(ajhg$id)))
tmp=fisher.test(matrix(c(a,c,117012-a,d),nrow = 2,byrow = T))
result[2,]=c("117K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(filea1$id),as.character(ajhg$id)))
tmp=fisher.test(matrix(c(a,c,61725-a,d),nrow = 2,byrow = T))
result[3,]=c("61K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file$id),as.character(ajhg$id)))
tmp=fisher.test(matrix(c(a,c,53425-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file1$id),as.character(ajhg$id)))
tmp=fisher.test(matrix(c(a,c,8544-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

write.csv(result,"20201203/AJHG.statis.csv",quote=F,row.names = F)

###220K ASM
c=length(intersect(as.character(con$id),as.character(asmid)))
d=dim(con)[1]-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",c,NA,NA)

a=length(intersect(as.character(filea$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,117012-a,d),nrow = 2,byrow = T))
result[2,]=c("117K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(filea1$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,61725-a,d),nrow = 2,byrow = T))
result[3,]=c("61K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,53425-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file1$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,8544-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file2$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,807-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file3$id),as.character(asmid)))
tmp=fisher.test(matrix(c(a,c,200-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"20201203/220K.ASM.statis.csv",quote=F,row.names = F)



###imprinted locis
library(ChIPpeakAnno)
imprint=read.table("E:/0 公共数据库差异情况/db_for_5hmc/imprint.hg19verion.bed",head=F,sep="\t")
names(imprint)=c("chr","start","end")
bed.imprint=toGRanges(imprint,format="BED",header=TRUE)

same.direction=read.csv("./20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",header=T)
same.direction=same.direction[same.direction$upup.num>1|same.direction$downdown.num>1,]
bed.same=data.frame(chr=same.direction$Chr,start=as.numeric(same.direction$Start)-1,end=as.numeric(same.direction$Start))
bed.same=toGRanges(bed.same,format="BED",header=TRUE)

same.direction1=read.csv("./20201124/same.direction.dowdown.upup.modify.csv",header=T)
same.direction1=same.direction1[same.direction1$group.DC.num>1,]
sel=which(same.direction1$pattern.no.dupl.DC=="up:down"|same.direction1$pattern.no.dupl.DC=="down:up")
same.direction1=same.direction1[-sel,]
same.direction1=tidyr::separate(same.direction1,unitID,into=c("Chr","Pos","Start","End"),sep=":")
same.direction1$id=paste(same.direction1$Chr,same.direction1$Pos,sep=":")
bed.same.direction1=data.frame(chr=same.direction1$Chr,start=as.numeric(same.direction1$Pos)-1,end=as.numeric(same.direction1$Pos))
bed.same.direction1=toGRanges(bed.same.direction1,format="BED",header=TRUE)

bed.filea=data.frame(chr=filea$Chr,start=as.numeric(filea$Start)-1,end=as.numeric(filea$Start))
bed.filea=toGRanges(bed.filea,format="BED",header=TRUE)

bed.filea1=data.frame(chr=filea1$Chr,start=as.numeric(filea1$Start)-1,end=as.numeric(filea1$Start))
bed.filea1=toGRanges(bed.filea1,format="BED",header=TRUE)

bed.file=data.frame(chr=file$Chr,start=as.numeric(file$Start)-1,end=as.numeric(file$Start))
bed.file=toGRanges(bed.file,format="BED",header=TRUE)

bed.file1=data.frame(chr=file1$Chr,start=as.numeric(file1$Start)-1,end=as.numeric(file1$Start))
bed.file1=toGRanges(bed.file1,format="BED",header=TRUE)

bed.file2=data.frame(chr=file2$Chr,start=as.numeric(file2$Start)-1,end=as.numeric(file2$Start))
bed.file2=toGRanges(bed.file2,format="BED",header=TRUE)

bed.file3=data.frame(chr=file3$Chr,start=as.numeric(file3$Start)-1,end=as.numeric(file3$Start))
bed.file3=toGRanges(bed.file3,format="BED",header=TRUE)

bed.con=data.frame(chr=con$Chr,start=as.numeric(con$Pos)-1,end=as.numeric(con$Pos))
bed.con=toGRanges(bed.con,format="BED",header=TRUE)

ol1=findOverlapsOfPeaks(bed.con,bed.imprint)
tcon=as.data.frame(ol1$peaklist$`bed.con///bed.imprint`)

ol2=findOverlapsOfPeaks(bed.filea,bed.imprint)
t2=as.data.frame(ol2$peaklist$`bed.filea///bed.imprint`)

ol3=findOverlapsOfPeaks(bed.filea1,bed.imprint)
t3=as.data.frame(ol3$peaklist$`bed.filea1///bed.imprint`)

ol4=findOverlapsOfPeaks(bed.file,bed.imprint)
t4=as.data.frame(ol4$peaklist$`bed.file///bed.imprint`)

ol5=findOverlapsOfPeaks(bed.file1,bed.imprint)
t5=as.data.frame(ol5$peaklist$`bed.file1///bed.imprint`)

ol6=findOverlapsOfPeaks(bed.file2,bed.imprint)
t6=as.data.frame(ol6$peaklist$`bed.file2///bed.imprint`)

ol7=findOverlapsOfPeaks(bed.same,bed.imprint)
t7=as.data.frame(ol7$peaklist$`bed.same///bed.imprint`)

ol8=findOverlapsOfPeaks(bed.file3,bed.imprint)
t8=as.data.frame(ol8$peaklist$`bed.file3///bed.imprint`)

ol9=findOverlapsOfPeaks(bed.same.direction1,bed.imprint)
t9=as.data.frame(ol9$peaklist$`bed.same.direction1///bed.imprint`)

c=dim(tcon)[1]
d=dim(con)[1]-c
result=data.frame(matrix(,nrow = 9,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",dim(tcon)[1],NA,NA)

tmp=fisher.test(matrix(c(dim(t2)[1],c,117012-dim(t2)[1],d),nrow = 2,byrow = T))
result[2,]=c("117K",dim(t2)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t3)[1],c,61725-dim(t3)[1],d),nrow = 2,byrow = T))
result[3,]=c("61K",dim(t3)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t4)[1],c,53425-dim(t4)[1],d),nrow = 2,byrow = T))
result[4,]=c("53425",dim(t4)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t5)[1],c,8544-dim(t5)[1],d),nrow = 2,byrow = T))
result[5,]=c("8544",dim(t5)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t6)[1],c,807-dim(t6)[1],d),nrow = 2,byrow = T))
result[6,]=c("807",dim(t6)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t7)[1],c,3969-dim(t7)[1],d),nrow = 2,byrow = T))
result[7,]=c("3969",dim(t7)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t8)[1],c,200-dim(t8)[1],d),nrow = 2,byrow = T))
result[8,]=c("200",dim(t8)[1],tmp$estimate,tmp$p.value)

tmp=fisher.test(matrix(c(dim(t9)[1],c,1195-dim(t9)[1],d),nrow = 2,byrow = T))
result[9,]=c("1195",dim(t9)[1],tmp$estimate,tmp$p.value)

write.csv(result,"20201203/imprinted.locis.statis.csv",quote=F,row.names = F)